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The  internal  waves  in  the  south  of  the  Strait  of  Messina  (Italy)  are  studied  using  observational  data  and 
numerical  simulations.  The  observational  data  consisted  of  SAR  images,  XBT  probes,  CTD  yoyos  and 
thermistor  string  profiles  from  the  Coastal  Ocean  Acoustic  Changes  at  High  Frequencies  (COACH06)  cruise. 
The  numerical  Lamb  model  was  used  to  solve  the  fully  nonlinear,  nonhydrostatic  Boussinesq  equations  on  an 
/-plane.  The  model  is  2.5  dimensional  with  spatial  variation  in  a  vertical  plane  extending  along  depth  and  in 
along-stream  direction,  and  neglects  derivatives  in  cross-stream  direction. 

Eleven  out  of  fifteen  SAR  images  contained  internal  wave  events  for  the  month  of  October  2006.  From  these 
images,  we  estimated  some  of  the  internal  wave  characteristics:  a  least  square  fit  of  the  front  positions  of  the 
southwards  propagating  internal  wave  trains  gives  a  propagation  speed  of  1.00  m  s_  1  and  a  time  of  release  of 
the  bore  of  about  2  h  after  maximum  northwards  tidal  flow  at  Punta  Pezzo.  Nonhydrostatic  dispersion  leads 
to  a  wavelength  increase  between  the  first  two  internal  solitary  waves  of  the  trains  of  40  m  km- 1. 

Our  in-situ  data  were  used  to  initialize  and  evaluate  the  Lamb  model.  An  EOF  analysis  was  applied  to  the  data 
and  the  model  outputs.  The  first  three  empirical  functions  contain  over  99%  of  the  variability.  The  data  and 
the  model  results  are  in  very  good  qualitative  and  quantitative  agreement,  giving  a  propagation  depth  of  the 
internal  solitary  wave  train  around  90  m  with  a  pycnocline  oscillating  from  80  to  130  m. 

Using  the  first  two  EOFs,  an  original  method  for  detecting  internal  waves  was  developed  by  analyzing  the 
scatter  diagram  of  the  first  EOF  projection  coefficient  versus  the  second  EOF  projection  coefficient.  This 
distribution  has  a  specific  “crescent  shape”  showing  that  the  first  two  projection  coefficients  are  not 
independent  in  the  presence  of  internal  wave  events.  The  “crescent  shape”  signature  of  the  internal  waves 
can  be  used  as  an  internal  wave  detector. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  Strait  of  Messina  is  a  rich  dynamic  area  where  many  ocean 
features  can  be  encountered.  It  separates  the  Italian  Peninsula  from 
the  island  of  Sicily  and  is  a  natural  connection  between  the  Tyrrhenian 
Sea  in  the  north  and  the  Ionian  Sea  in  the  south  (Fig.  1 ).  The  strait  is  a 
narrow  channel  whose  smallest  cross-sectional  area  is  0.3  km2  in  the 
sill  region.  The  sill  raises  to  within  80  m  of  the  surface  at  the 
shallowest  point.  The  water  depth  increases  gently  in  the  northern 
part  of  the  strait  whereas  the  slope  is  steeper  in  the  south  basin 
(15  km  away  from  the  sill,  the  depth  is  400  m  in  the  north  basin,  and 
800  m  in  the  south).  Throughout  the  year,  two  water  masses  are 
encountered  in  the  strait:  the  Tyrrhenian  surface  water  and  the  colder 
and  saltier  Levantine  Intermediate  Water  from  the  Ionian  Sea. 

Although  the  Strait  of  Messina  has  been  known  since  ancient  times 
as  an  area  of  strong  currents  and  vortices  (represented  by  the  two 
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mythological  monsters  Scylla  and  Charybdis  in  Homer  Odyssey, 
800  B.C.),  more  than  2000  years  elapsed  before  modern  scientific 
oceanographic  measurements  were  carried  out  in  this  area.  The  French 
vice-consul  Ribaud  gave  a  fairly  detailed  description  of  the  currents  in 
the  strait  in  1824,  but  it  was  not  until  1922-23  that  an  extensive 
oceanographic  survey  was  made  by  the  Italian  oceanographer  Vercelli 
(1925)  with  the  research  vessel  Marsigli.  The  data  collected  are  still 
considered  to  be  the  most  detailed  and  systematic  data  from  the  Strait 
of  Messina  and  were  used  and  processed  by  Hopkins  et  al.  in  1984. 
What  really  initialized  the  intense  study  of  the  Strait  of  Messina  was  the 
first  satellite  observation  of  internal  waves  made  by  SEASAT  on  15 
September  1978  (Alpers  and  Salusti,  1983).  In  the  following  years, 
oceanographic  campaigns  were  carried  out  to  measure  internal  solitary 
waves  north  and  south  of  the  strait.  In  the  following  years,  a 
considerable  number  of  scientific  studies  of  internal  solitary  waves  in 
the  Strait  of  Messina  were  made  through  in-situ  and  satellite 
observations  (Alpers  and  Salusti,  1983;  Hopkins  et  al.,  1984;  Griffa 
et  al.,  1986;  Marullo  and  Santoleri,  1986;  Sapia  and  Salusti,  1987; 
Gerkema,  1994;  Brandt  et  al.,  1997;  Brandt  et  al.,  1998). 

From  the  measurements  carried  out,  a  current  distribution  has 
emerged.  The  dominating  role  of  semidiurnal  tides  was  described  in 
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Fig.  1.  Bathymetry  of  the  Strait  of  Messina.  The  points  show  the  measurement  positions:  blue  dot  is  the  thermistor  string  position,  pink  dots  indicate  the  CTD  stations,  and  black  dots 
represent  XBT  probes. 


detail,  and  the  influence  of  baroclinicity  on  the  formation  of  the  velocity 
field  in  the  strait  was  qualitatively  determined  in  Hopkins  et  al.  (1984) 
and  Mosetti  (1988).  At  the  sill,  a  stationary  surface  current  is  flowing  in 
the  southwards  direction  («0.10  m  s-1  time-averaged  velocity)  and  a 
bottom  current  is  flowing  in  the  opposite  direction  towards  the 
Tyrrhenian  Sea  (^0.13  m  s-1  time-averaged  velocity).  The  Ionian 
water  coming  from  the  south  is  saltier  and  heavier,  the  relative  density 
difference  is  of  the  order  of  0.001.  Superimposed  on  these  stationary 
currents  are  the  tidal  currents  originating  from  the  co-oscillation  of  the 
water  masses  of  the  strait  with  the  tides  of  the  adjacent  seas.  During 
maximum  tidal  flow  from  the  Tyrrhenian  Sea  to  the  Ionian  Sea  (rema 
scendente)  and  inversely  (rema  montante),  only  one  type  of  water  is 
encountered  over  the  sill,  thus  annihilating  the  two-layer  structure  over 
the  sill:  the  tidal  barotropic  aspect  of  water  motion  was  analyzed  by 
Defant  (1961).  The  time-averaged  interface  between  these  two  water 
masses  is  at  a  depth  of  150  m  in  the  adjacent  basins  but  at  30  m  around 
and  above  the  strait  which  means  an  uplift  of  120  m  from  the  adjoining 
basins  (Bignami  and  Salusti,  1990).  Although  tidal  sea  surface  displace¬ 
ments  are  very  small  in  the  strait  of  Messina  (the  order  of  10  cm),  Defant 
( 1940)  explained  the  existence  of  strong  currents  in  the  strait  by  the  fact 
that  the  tides  for  the  open  strait  boundaries  oscillate  almost  in 
opposition,  thus  leading  to  a  sea  surface  slope  of  1-2  cm  km-1 
generating  strong  currents  (Mosetti,  1988)  as  high  as  3  m  s_  1  in  the  sill 
region  during  spring  tides  (Bignami  and  Salusti,  1990). 

The  strong  barotropic  tidal  flow  over  a  steep  bathymetry  in  a  stably 
stratified  environment  represents  the  ideal  recipe  for  internal  wave 
generation  (Zeilon,  1912).  Baines  (1973, 1982)  has  demonstrated  the 
generation  process  analytically.  Hydrodynamic  modelling  shows  that 
tidal  flow  over  steep  topography  creates  an  interfacial  depression  of 
the  pycnocline  and  generates  an  internal  bore  (Del  Ricco,  1982;  Lamb, 
1994;  Mosetti,  1988;  Brandt  et  al.,  1997;  Warn  Varnas  et  al.,  2003, 
2005, 2007).  In  the  Strait  of  Messina,  during  rema  montante  when  the 
heavier  Levantine  Intermediate  Water  crosses  the  sill,  the  pycnocline 
moves  upwards  at  the  sill  and  is  depressed  in  the  northern  basin.  The 
depression  generates  a  southwards  and  a  northwards  propagating 


bore.  The  northwards  propagating  bore  steepens  and  disintegrates 
into  solitary  waves  whereas  the  southwards  propagating  bore  is 
stopped  by  the  sill.  When  the  semidiurnal  tide  reverses  to  rema 
scendente,  the  southwards  propagating  internal  depression  under¬ 
goes  a  hydraulic  jump  over  the  sill  and  into  the  southern  basin  where 
it  propagates  away  from  the  sill.  As  it  propagates,  nonlinear  effects 
steepen  its  leading  edge  until  it  disintegrates  into  a  series  or  “train”  of 
interfacial  nonlinear  short  internal  waves.  This  disintegration  is  due  to 
effects  of  frequency  and  amplitude  dispersions,  as  well  described  in 
Warn  Varnas  et  al.  (2007).  Internal  waves  in  the  Strait  of  Messina  are 
characterized  by  a  propagation  speed  between  0.80  to  1.00  m  s_1, 
with  oscillations  of  temperature  up  to  2  °C  amplitude.  When  the  train 
is  well  formed,  4  to  10  internal  solitary  waves  per  train  can  be 
observed  with  periods  ranging  from  8  to  30  min,  covering  an  average 
total  duration  of  about  2  h  (Alpers  and  Salusti,  1983). 

In  this  paper,  we  present  an  original  way  of  characterizing  the 
internal  waves  in  the  Strait  of  Messina  by  using  their  multimodal 
structure.  For  October  2006,  fifteen  SAR  images  gave  the  internal 
waves  characteristics  which  are  used  to  tune  a  Lamb  numerical  model 
for  the  Strait  of  Messina.  The  model  was  initialized  from  in-situ  data 
and  numerical  outputs  are  compared  with  the  COACH06  cruise 
oceanographic  measurements  (late  October-beginning  of  November 
2006)  by  applying  an  Empirical  Orthogonal  Function  (EOF)  analysis. 
Using  the  crescent-shape  scatter  plot  formed  by  the  first  and  second 
EOF  projection  coefficients,  it  was  possible  to  characterize  the  presence 
of  internal  waves  south  of  the  Strait.  Section  2  describes  the  dataset  and 
the  EOF  analysis  of  the  in-situ  data.  Section  3  presents  the  model  and 
Section  4  discusses  the  scatter  diagram  results.  Section  5  contains  a 
summary  and  conclusions. 

2.  Dataset  exploitation 

The  dataset  consisted  of  SAR  images  (October  2006)  and  in-situ 
hydrographical  data  (4-11  November  2006)  such  as  XBT  probes,  CTD 
yoyos  and  a  thermistor  chain. 
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Internal  Waves 
in  the  Strait  of  Messina 
21  October  2006  16:48  GMT 

1 5g25‘E  15°30'E  15°35'E  15°40‘E  15°45'E 


Fig.  2.  SAR  image  from  RADARSAT1  satellite  with  a  clear  internal  solitary  waves  train  signature  seen  on  21  October  2006, 16:48  UTC. 


2  A.  SAR  images 

The  propagating  solitary  waves  are  seen  as  isopycnal  depressions 
(Osborne  and  Burch,  1980).  When  such  interfacial  depressions 
propagate,  they  generate  a  downstream  flow  convergence  that 
increases  the  surface  roughness  and  upstream  flow  divergence  that 
tends  to  smooth  the  surface  roughness.  SAR  are  most  sensible  to  such 
surface  roughness  and  are  able  to  picture  internal  waves  as  trains  of 
bright  (convergence)  and  dark  (divergence)  streaks  (Alpers  and 
Salusti,  1983;  Apel  and  Jackson,  2004).  The  first  time  that  tidal  bores 
have  been  identified  on  a  SAR  image  in  the  Strait  of  Messina  was  on  15 
September  1978  (Alpers  and  Salusti,  1983).  Since  then,  the  SAR  data 
have  been  widely  used  in  internal  wave  studies  in  the  Strait  of  Messina 
(Alpers  et  al.,  1996;  Brandt  et  al.,  1997).  Global  Ocean  Associates 
(GOA),  founded  by  Dr.  John  R.  Apel  in  1996  and  funded  by  the  Office  of 
Naval  Research,  made  a  worldwide  Atlas  of  Internal  Solitary-like 
Waves  (GOA,  2004)  using  SAR  images.  In  this  Atlas,  internal  solitary 
wave  characteristics  were  computed  for  the  Strait  of  Messina  and  the 
adjacent  sea  areas,  by  using  SAR  images  from  ERS-1/2  satellites  for  a 
period  ranging  from  beginning  of  1991  to  December  1995.  On  the 
whole  period,  77  manifestations  of  internal  waves  could  be  delineated 
out  of  160  overflights.  That  study  showed  that  sea  surface  manifesta¬ 
tions  of  internal  solitary  wave  trains  are  observed  more  frequently 


during  summer  when  a  strong  seasonal  pycnocline  is  known  to  be 
present.  Furthermore,  sea  surface  manifestations  of  southwards 
propagating  internal  solitary  wave  trains  are  more  frequent  than 
northwards  propagating  ones.  For  the  month  of  October,  on  average, 


Table  1 

SAR  image  characteristics. 


Date 

Time  (UTC) 

At  (IM/MNF) 

d.Sill  (m) 

n 

Length  (m) 

d.1/2  (m) 

3rd  Oct  2006 

9:02 

08:17 

22703 

5 

3784 

946 

5th  Oct  2006 

5:01 

15:25 

46637 

7 

6168 

1852 

6th  Oct  2006 

9:08 

18:08 

59757 

5 

6000 

2571 

8th  Oct  2006 

21:01 

17:06 

54769 

12 

17700 

2187 

11th  Oct  2006 

16:43 

10:54 

31634 

2 

1500 

1416 

15th  Oct  2006 

20:41 

09:52 

29142 

3 

3429 

1714 

21st  Oct  2006 

16:48 

13:58 

41651 

12 

10769 

2051 

21st  Oct  2006 

20:52 

18:02 

56456 

10 

13021 

2532 

6:00 

14142 

1 

X 

X 

24th  Oct  2006 

17:00 

12:52 

35833 

4 

3789 

1428 

24th  Oct  2006 

20:58 

16:50 

51210 

2 

2032 

2034 

‘Date’  and  ‘time’  of  the  SAR  image.  ‘At  (IM/MNF)’  is  the  time  difference  between  the  SAR 
image  and  the  previous  occurrence  of  maximum  northwards  flow  in  P.  Pezzo,  ‘d.Sill’  is 
the  distance  between  the  observed  internal  wave  train  and  the  sill,  ‘n’  is  the  number  of 
solitary  waves  in  the  train,  ‘length’  is  the  length  of  the  whole  train,  ‘d.1/2’  is  the  spatial 
separation  between  the  first  and  second  waves. 
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there  are  10  observations  of  internal  solitary  wave  trains  out  of  14 
overflights  south  of  the  Strait  of  Messina. 

In  our  dataset,  15  SAR  images  for  October  2006  were  obtained  from 
the  ENVISAT  observation  satellite  (launched  in  2002  by  the  European 
Space  Agency,  ESA)  and  from  the  RADARSAT-I  satellite  (launched  in 
1995  by  the  Canadian  Space  Agency,  CSA).  Eleven  of  them  presented 
internal  solitary  wave  train  signatures  (Fig.  2).  Table  1  summarizes,  for 
each  SAR  image  showing  an  internal  wave  event,  the  date  and  time  of 
the  image,  the  time  difference  between  the  SAR  image  and  the 
previous  occurrence  of  maximum  northwards  flow  in  P.  Pezzo,  the 
distance  of  the  first  internal  solitary  wave  from  the  sill,  the  number  of 
solitary  waves  in  the  train,  the  length  of  the  whole  train  and  the 
spatial  separation  between  the  first  and  second  waves.  The  internal 
wave  train  propagation  velocity,  the  time  of  release  of  the  bore  and  the 
spatial  separation  between  the  first  two  internal  waves  of  the  trains 
have  been  computed  from  these  data  in  Sections  2.1.1  and  2.1.3.  These 
parameters  give  a  statistical  representation  of  the  distinctive  south¬ 
wards  propagating  internal  waves  formed  at  the  sill. 


2.1.1.  Internal  wave  propagation  velocity  and  time  of  release  of  the  bore 
A  time  versus  distance  diagram  delineating  the  propagation  of 
internal  bores  is  shown  in  order  to  estimate  the  propagation  velocity  of 
the  internal  solitary  wave  train  and  the  time  of  release  of  internal  bores 
from  the  sill  (Fig.  3  top).  The  “distance”  is  measured  along  a  mean  path 
between  the  sill  and  the  centre  of  the  wave  fronts,  which  is  inferred 
from  the  available  SAR  images;  “time”  refers  to  the  time  between  the 
occurrence  of  the  previous  maximum  northwards  flow  in  P.  Pezzo  and 
the  moment  when  SAR  image  was  taken.  The  line  is  a  least  square 
fit  that  yields  an  estimated  propagation  speed  of  1.00  m  s_  l  The  time 
of  release  of  the  bore  at  the  sill  is  about  2  h,  which  is  in  excellent 
agreement  with  the  Atlas  of  Internal  Waves  (GOA,  2004),  which  with 
the  same  methodology  yields  a  propagation  velocity  of  0.91  m  s_  1  and 
a  time  of  release  of  1  to  5  h.  We  have  described  in  the  introduction  that 
the  bore  was  generated  by  an  interfacial  depression  at  the  lee  side  of 
the  sill  by  the  northwards  tidal  flow,  and  is  released  into  the  southern 
basin  when  the  tidal  flow  weakens  and  reverses  southwards.  The 
time  of  release  of  the  bore  is  then  the  time  between  the  maximum 


Time  after  max  North  flow  in  P. Pezzo  versus  IW  position 


First  wavelength  versus  IW  position 


South  (Ionian  Sea)  <-  Distance  from  sill  (km)  ->  Sill 

Fig.  3.  Study  of  the  11  SAR  images  showing  internal  wave  trains  south  of  Messina  Strait  (Ionian  Sea).  Top,  space-time  diagram  showing  the  propagation  of  the  front  of  internal  waves. 
Middle,  number  of  internal  solitary  waves  per  train  versus  the  IW  train  position  in  the  south  of  the  strait.  Bottom,  spatial  separation  between  the  first  two  internal  solitary  waves  of 
each  train  as  a  function  of  the  IW  train  position.  “Distance  from  sill”  (x-axis)  is  the  IW  train  position  represented  as  the  distance  between  the  first  wave  and  the  sill.  South  is 
considered  to  be  on  the  left  side,  north  on  the  right.  “Time”  is  the  temporal  separation  between  the  maximum  northwards  flow  in  Punta  Pezzo  and  the  time  of  the  SAR  image.  The 
lines  are  least  square  fits. 
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northwards  flow  in  P.  Pezzo  and  the  slack  water  corresponding  to  the 
tidal  flow  reversal  when  the  bore  undergoes  its  hydraulic  jump  over 
the  sill.  Looking  at  tide  tables  for  the  month  of  October  2006,  this  time 
of  release  varies  from  2.5  to  4  h,  fitting  in  the  interval  of  1  to  5  h 
mentioned  in  the  Atlas  of  Internal  Waves  (2004). 

2.1.2.  Number  of  waves  versus  the  distance  away  from  the  sill 

The  second  graph  of  Fig.  3  (middle)  illustrates  the  disintegration  of 
the  internal  bore  into  a  train  of  internal  solitary  waves.  It  represents 
the  number  of  internal  solitary  waves  per  train  versus  the  distance  to 
the  strait.  Only  one  SAR  image  shows  an  event  close  to  the  strait,  that 
is  one  observation  at  14  km  from  the  sill.  This  observation  can  be 
interpreted  as  a  bore  that  has  not  yet  disintegrated.  It  is  only  at  least 
23  km  away  from  the  sill  that  the  number  of  waves  increases,  showing 
when  the  bore  disintegrates  into  internal  solitary  waves.  This  means 
that  the  bore  has  to  travel  a  very  long  distance  before  disintegrating 
into  a  wave  train.  This  effect  is  mainly  due  to  the  strong  increase  in  the 
bottom  depth  south  of  the  sill  of  the  strait  as  discussed  in  Griffa  et  al. 
(1986). 

2.1.3.  Spatial  separation  between  the  first  two  waves  for  a  well  formed 
internal  wave  train 

Fig.  3  (bottom)  shows  the  spatial  separation  between  the  first  two 
internal  waves  of  the  propagating  wave  trains  as  a  function  of  the 


distance  from  the  sill.  The  first  wave  of  the  train  has  bigger  amplitude 
and  velocity  than  the  second  wave,  meaning  that  the  wavelength 
between  the  first  two  internal  solitary  waves  increases  as  the  train 
propagates.  In  the  Atlas  of  Internal  Waves  (2004),  the  spatial 
separation  ranges  from  500  m  to  1900  m.  For  our  set  of  images,  the 
wavelength  ranges  from  1300  m  (30  km  south  of  the  sill)  to  2500  m 
(60  km  south  of  the  sill),  meaning  that  the  spatial  separation  between 
the  first  two  waves  increases  by  40  m  km-1  as  the  train  travels 
southwards.  Still  this  figure  must  be  considered  cautiously  as  it  is  very 
difficult  to  measure  exactly  the  spatial  separation  of  the  first  two 
waves  because  of  the  imperfect  circular  shape  of  the  wave  fronts  (see 
wave  front  on  Fig.  2). 

2.2.  In-situ  data 

During  the  COACH06  cruise  (4-11  November  2006),  the  French 
R/V  Beautemps-Beaupre  made  CTD  and  XBT  measurements.  One 
thermistor  chain  was  also  deployed  (Fig.  1). 

49  XBT  probes  were  launched  between  4  and  7  November  in  the 
southern  basin.  These  XBT  measurements  provided  temperature 
profiles,  and  density  was  computed  from  Unesco  equation  of  state. 
The  density  profiles  show  great  variability,  illustrating  the  presence  of 
internal  solitary  waves  (Fig.  4  top).  The  XBT  measurements  were  used 
to  check  whether  the  distance  from  the  coasts  would  influence  the 


XBT  density  profiles 


Density 

CTD2  Density  profiles 


Density 

Fig.  4.  Top,  density  profiles  from  all  XBT  probes  (see  Fig.  1  for  locations)  in  the  south  of  the  strait  (Ionian  Sea)  illustrating  the  variability  in  the  water  column  due  to  the  impact  of 
internal  waves.  Bottom,  density  profiles  measured  by  the  lowered  CTD  probe  at  location  (15°34/04"  E,  38°03'04"  N).  The  pycnocline  depth  has  a  25  m  amplitude  variability  illustrating 
the  occurrence  of  an  internal  wave  event.  Density  units  are  “sigma-t”. 
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Interpolated  density  from  thermistor  string  measurement 


5-Nov-06  (1200UTC)  6-Nov-06  (1200UTC)  7-Nov-06 


Fig.  5.  Density  computed  from  thermistor  string  measurements  from  4  November  2006  16:48  UTC  to  7  November  2006  07:12  UTC.  The  28.5  density  contour  (corresponding  to  a 
temperature  of  16  °C)  is  plotted  to  underline  the  temperature  variability  due  to  occurrence  of  internal  waves.  Density  units  are  “sigma-t”. 


depth  of  the  mixed  layer,  since  one  would  expect  more  mixing  and 
more  meteorological  impact  closer  to  the  coast.  So,  the  measurements 
were  divided  in  three  zones  (west,  middle,  east  of  the  strait).  An 
average  profile  was  given  for  each  zone  (not  shown).  All  this  shows 
that  the  mixing  is  constant  from  the  Calabrian  coast  (east)  to  the 
middle  of  the  Strait  and  increases  closer  to  the  Sicilian  coast  (west). 
This  could  be  due  to  coastal  current  dynamics  on  the  eastern  coast  of 
Sicily  as  described  in  Bohm  et  al.  (1987).  The  tidal  mixing  in  Messina 
Strait  generates  a  water  mass  that  flows  up  to  100  km  southwards 
along  the  Sicilian  coast  and  up  to  10  km  offshore.  This  current  can 
become  unstable  and  lead  to  more  mixing.  On  the  SAR  image 
displayed  in  Fig.  2,  it  is  clearly  visible  that  the  sea  surface 
manifestations  are  stronger  near  the  Sicilian  coast  than  further 
offshore  or  close  to  the  Calabrian  coast.  The  amplitude  of  the  internal 
waves  is  larger  in  this  shallow  water  region  of  the  Sicilian  coast,  where 
the  internal  waves  are  topographically  guided  (Alpers  et  al.,  1996).  As 
the  measurements  were  done  close  to  the  Calabrian  coast  and  the 
model  runs  along  the  median  of  the  strait,  this  difference  is  neglected 
in  the  following. 

Lowered  CTD  probes  measurements  (yoyos)  were  done  at  three 
locations:  Station  1  (15°30'19"  E,  38°03'54"  N)  from  6  November 
08:12  UTC  to  6  November  12:50  UTC,  Station  2  (15°34/04"  E,  38°03' 
04"  N)  from  6  November  20:38  UTC  to  7  November  01 :47  UTC,  Station 
3  (15°30'19"  E,  38°03'01"  N)  from  7  November  13:49  UTC  to  7 
November  15:07  UTC.  Fig.  4  (bottom)  illustrates  the  density  profiles 
for  the  CTD  measurements  at  station  2.  The  low  temporal  sample  rate 
of  the  lowered  CTD  (one  profile  every  30  min)  makes  it  impossible  to 
detect  the  propagation  of  a  solitary  wave  train  as  the  lifting  and  return 
to  equilibrium  of  the  isopycnals  is  much  faster,  with  internal  solitary 
wave  periods  ranging  from  8  to  30  min  (Alpers  and  Salusti,  1983). 
However,  the  pycnocline  depth  varies  by  25  m  in  5  h  of  measurement 
clearly  indicating  some  internal  wave  feature. 

Finally,  as  done  and  validated  in  Sellschopp  ( 1997)  for  high  temporal 
resolution  observation  of  isothermal  displacements,  one  thermistor 
chain  was  bottom-moored  about  27  km  south  of  the  sill  of  strait,  at 
location  (15°38'57"  E,  38°02'24"  N),  from  4  November  16:48  UTC  to  the 
7  November  07:12  UTC,  with  a  one  minute  temporal  resolution  which 
allows  clear  observations  of  internal  solitary  wave  events.  The  chain 
consisted  of  10  sensors  attached  from  14  to  128  m  depth,  providing  a 
vertical  resolution  of  about  10  m.  Fig.  5  shows  the  density  computed 
from  the  temperature  measurement  of  the  thermistor  string  using  the 
Unesco  equation  of  state  (salinity  was  considered  constant  with  a  value 
of  38.10  psu).  The  pycnocline  depth  stands  around  95  m  for  the  first  day 
and  half,  and  raises  slowly  to  a  60  m  depth  after  6  November.  Internal 
solitary  wave  train  events  can  be  identified  every  half  day,  correspond¬ 


ing  to  the  tidal  period.  Table  2  summarizes  precisely  the  occurrence  of 
the  internal  wave  trains,  the  reversal  of  the  tide  from  rema  montante 
(northwards)  to  rema  scendente  (southwards).  On  average,  the  internal 
wave  trains  are  visible  6  h  after  the  tide  reversal  from  North  to  South. 
Considering  the  hypothesis  made  by  Alpers  and  Salusti  (1983)  that,  for 
southwards  propagating  internal  solitary  waves,  the  internal  bores  are 
generated  at  slack  waters  when  the  flow  reverses  from  north  to  south, 
the  average  speed  of  the  internal  waves  train  is  1.25  m  s_1,  which  is 
faster  than  the  propagation  speed  computed  from  the  SAR  images.  A 
difference  in  the  stratification,  in  the  mean  current  and/or  the  tidal 
cycles  between  the  months  of  October  and  November  could  explain  this 
difference.  Unfortunately  no  SAR  images  were  available  for  beginning  of 
November.  The  amplitudes  of  the  internal  solitary  waves  range  between 
40  and  50  m.  The  trains  have  two  to  four  internal  solitary  waves.  The  first 
and  second  waves  are  very  neat  in  the  data.  This  is  in  good  agreement 
with  the  analysis  of  the  SAR  images:  on  Fig.  3  (middle  graph),  the  least 
square  fit  gives  a  value  close  to  four  internal  waves  per  train  at  27  km 
south  of  the  sill,  at  the  thermistor  string  position. 

This  data  analysis  was  used  to  tune  the  Lamb  model  (third  section) 
to  infer  the  right  internal  wave  characteristics  by  focusing  on 
propagation  velocity,  number  of  internal  waves  per  train,  increase  in 
the  spatial  separation  of  the  waves  (as  obtained  from  the  analysis  of 
the  SAR  images),  and  frequency  and  amplitude  of  the  internal  solitary 
waves  (as  observed  from  in-situ  data). 


Table  2 

Internal  wave  trains  (IWT)  identified  on  the  thermistor  string  measurement. 


Internal  wave  IWT  1 
trains 

IWT  2 

IWT  3 

IWT  4 

IW  5 

Occurrence 

Nov.  4th  2006  Nov.  5th  2006  Nov.  5th  2006  Nov.  6th  2006  Nov.  7th  2006 

of  IWT 

23:31  UTC 

11:55  UTC 

23:40  UTC 

13:15  UTC 

00:15  UTC 

Previous  tide 

Nov.  4th  2006  Nov.  5th  2006 

Nov  5th  2006 

Nov.  6th  2006 

Nov.  6th  2006 

reversal 

17:02  UTC 

06:05  UTC 

17:39  UTC 

06:48  UTC 

18:14  UTC 

(slack 

waters) 

At  (REV-IWT) 

6:29 

5:50 

6:01 

6:27 

6:01 

As  the  bore  generated  at  “rema  montante”  on  the  northern  side  of  the  strait  undergoes  a 
hydraulic  jump  over  the  sill  to  the  southern  basin  when  the  tide  reverses  southwards 
(rema  scendente),  the  time  difference  between  the  SAR  image  and  the  tidal  reversal  is 
computed  in  order  to  get  the  propagation  speed  of  the  internal  wave  train  in  the  south 
basin.  The  first  line  indicates  the  date  and  time  of  the  IWT  event,  the  second  line 
indicates  the  time  of  tidal  reversal  from  north  to  south  (slack  waters)  previous  to  the 
SAR  image.  The  last  line  shows  the  time  difference  between  the  previous  tide  reversal 
and  the  occurrence  of  the  IWT. 
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2.3.  Empirical  orthogonal  function  analysis  of  thermistor  chain  data 

An  EOF  analysis  of  the  time  series  of  density  obtained  from  in-situ 
data  and  numerical  simulations  was  carried  out  to  gain  more 
information  on  the  internal  wave  structure. 

The  empirical  orthogonal  function  (EOF)  decomposition  is  one  of  a 
larger  class  of  inverse  techniques  and  is  equivalent  to  a  data  reduction 
method  (Emery  and  Thomson,  1997).  It  determines  the  main  spatial 
patterns  of  variability  as  well  as  their  variation  in  time  and  gives  a 
measure  of  importance  of  each  pattern.  Each  EOF  mode  represents  a 
standing  oscillation  pattern  (Preisendorfer,  1988).  The  projection  of 
the  data  onto  the  EOF  modes  gives  an  estimate  of  how  the  pattern 
oscillates  in  time.  These  projections  are  called  principal  component 
time  series  or  expansion/projection  coefficients  of  the  EOFs.  Most  of 
the  variance  is  usually  in  the  first  few  orthogonal  functions  and  the 
EOFs  can  be  employed  as  a  filter  to  avoid  experimental  or  physical 
noises  and  consider  only  the  most  important  scales  of  variability. 

In  our  study,  the  EOF  analysis  was  applied  on  the  vertical  density 
anomaly  obtained  from  the  thermistor  string  measurements,  so  that 
the  density  o[z,  t)  could  be  written  as: 

M 

CT(Z,t)  =  (<7(Z))  +  “n(t)</>„(Z)  ("•) 

n  =  1 

where  <cr(z)>  is  the  time-averaged  vertical  density  profile,  M  is  the 
number  of  considered  EOFs,  an(t)  is  the  amplitude  (or  ‘expansion 
coefficient’)  of  the  nth  orthogonal  mode  at  time  t  and  c pn(z)  the  nth 
EOF.  This  equation  says  that  the  time  variation  of  the  dependent  scalar 
variable  a(z,t )  at  each  depth  results  from  the  linear  combination  of  M 
spatial  functions  4>n,  whose  amplitudes  are  weighted  by  M  time- 
dependent  coefficients  an(t).  The  weights  an(t)  tell  how  the  spatial 
modes  f>n  vary  with  time.  As  there  are  many  terminologies  for  EOF 
analysis,  we  insist  on  the  fact  that  the  use  of  “EOF”  refers  to  the  vertical 
“spatial”  pattern  4>n  and  that  “expansion  coefficient”  (or  principal 
component,  or  projection  coefficient)  refers  to  the  “temporal 
patterns”  or  weights  an(t). 


With  the  inherent  efficiency  of  this  statistical  description,  a  very 
few  empirical  modes  are  needed  to  describe  the  fundamental 
variability  in  a  very  large  data  set  and  so  only  three  modes  were 
considered  here.  The  first  two  EOFs  are  plotted  in  Fig.  6.  They  are  very 
clean  and  smooth,  representing  a  dynamical  pattern  consisting  of  a 
background  stratification  interacting  with  internal  wave  occurrence. 
Moreover,  most  of  the  variance  of  the  measurements  is  in  the  first  two 
EOFs  (95%  and  4.5%). 

The  amplitude  of  the  first  EOF  does  not  change  sign  over  the  water 
column  and  has  one  extremum.  Therefore,  all  the  isopycnals  oscillate 
in  the  same  direction  (Vazquez  et  al.,  2006)  with  maximum 
displacement  at  85-90  m.  It  represents  a  vertical  displacement  of 
the  pycnocline.  The  second  EOF  has  two  extrema,  as  a  consequence  of 
the  orthogonality  constraint,  and  changes  sign  close  to  the  first  EOF 
maximum  depth.  This  indicates  that  isopycnals  are  moving  in  opposite 
directions  above  and  below  the  pycnocline  (Vazquez  et  al.,  2006);  it 
can  be  identified  as  a  first  baroclinic  mode.  It  represents  a  change  in 
slope  of  the  density  gradient.  Fig.  7  shows  one  density  profile 
projected  on  the  first  EOF  mode  <(j(z)>  +  ah(tn)</>i  (left)  for  several 
random  values  of  the  time-dependent  coefficient  a i  and  the  same 
density  profile  projected  on  the  second  EOF  (cr{z))  +  a2{tn)(j) 2  (right) 
for  several  values  of  time  of  a2  to  illustrate  the  previous  description  of 
the  modes.  As  these  profiles  are  projected  on  a  single  EOF,  they  are  not 
all  realistic  (especially  when  projected  on  EOF2)  but  this  illustration 
helps  understanding  the  variability  of  the  first  two  modes.  The  depth 
of  the  pycnocline  increases/deepens  with  higher  positive  values  of 
whereas  the  density  gradient  becomes  steeper,  even  negative  with 
higher  values  of  a2.  Looking  at  Fig.  8,  we  see  that  the  values  of  the 
expansion  coefficients  are  clearly  related  to  internal  wave  occur¬ 
rences,  both  coefficients  assuming  much  larger  absolute  values  during 
internal  wave  events. 

Interpreting  EOFs  must  always  be  done  cautiously  as  empirical 
modes  do  not  correspond  necessarily  to  true  dynamical  characteristics 
or  modes  of  physical  behaviour  (Dommenget  and  Latif,  2002). 
Fortunately,  by  studying  the  data,  it  is  clearly  seen  that  our  computed 
EOF  do  agree  very  well  with  the  dynamical  modes  of  internal  waves. 


EOF1 


EOF2 


Amplitude  Amplitude 

Fig.  6.  Amplitude  of  first  two  EOFs  cf>i  and  cf>2  calculated  from  XBT  (black  dots),  thermistor  string  (grey  line)  and  model  (black  line)  densities.  The  model  fits  quantitatively  and 
qualitatively  well  the  EOFs  of  the  measurements.  Only  the  first  two  EOF  are  represented  as  they  give  over  99%  of  the  variance  (95%  and  4.5%). 
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Profile  projected  on  EOF1  only 
(alpha2=0) 


Density 


Profile  projected  on  EOF2  only 
(alpha1=0) 


Density 


Fig.  7.  On  the  left,  one  density  profile  is  projected  on  EOF1  <j>!  ( a2{t )  =  0)  for  several  random  values  of  On  the  right,  the  same  density  profile  is  projected  on  EOF2  cj>2  (aii(t)  =  0) 
for  several  random  values  of  a2(t).  The  values  of  each  considered  expansion  coefficient  are  displayed  in  the  legends.  The  first  EOF  mode  shows  a  vertical  translation  of  the  pycnocline 
whereas  the  second  EOF  mode  illustrates  the  change  in  slope  of  the  density  gradient. 


As  most  of  the  variance  was  contained  in  the  first  two  spatial  modes 
and  as  the  expansion  coefficients  have  a  similar  behaviour,  we  looked  at 
the  scatter  diagram  formed  by  plotting  the  first  expansion  coefficient  cti 


versus  the  second  expansion  coefficient  a2  in  order  to  see  if  they  could 
be  related  in  the  presence  of  an  internal  wave  event.  The  result  reveals 
a  crescent  shaped  distribution  (Fig.  9)  meaning  that  the  two  EOF 


EOF  expansion  coefficients  for  thermistor  chain  density 


Fig.  8.  Expansion  coefficients  for  the  first  two  EOF  modes.  The  first  expansion  coefficient  clearly  represents  the  vertical  oscillation  of  the  pycnocline  in  time  (compared  with  density 
contour  28.5  on  Fig.  5).  Both  coefficients  have  stronger  (absolute)  values  during  internal  wave  events. 
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Thermistor  data  scatter  diagram 
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Fig.  9.  Scatter  diagram  for  the  thermistor  string  data.  The  scatter  diagram  is  obtained  by  plotting  the  first  two  EOF  coefficients  a i  and  a2  versus  one  another.  The  “crescent  shape”  of 
the  scatter  diagram  shows  that  the  two  main  EOF  expansion  coefficients  are  not  independent  in  the  presence  of  internal  waves.  The  “middle”  of  the  cloud  of  points  represents  the 
interval  between  two  internal  wave  events  when  the  profile  is  close  to  its  mean  position.  The  asymptotic  behavior  of  the  cloud  of  points  for  high  values  of  the  expansion  coefficients  is 
representative  of  internal  wave  events.  The  black  dots  indicate  the  location  of  the  four  initial  density  profile  coefficients  for  chapter  4.2.  Profiles  A,  B  and  C  are  taken  out  of  the  cloud  of 
points  whereas  profile  D  is  taken  well  inside  the  cloud  of  points,  at  the  bottom  of  the  “crescent”. 


projection  coefficients  are  not  independent:  a2(t)  =/[aq(t)).  The 
density  can  then  be  fully  described  as: 

CT(z,t)  =  (a(z))  +  a,(t)<Mz)  +  a2(t)<l>2(z),  (2) 

=  Mz)>  +  c*i  fa  +f(ax)<t>2-  (3) 

When  looking  at  the  cloud  of  points  on  Fig.  9,  two  regimes  can  be 
distinguished.  In  the  centre  of  the  distribution  (|ah|<  0.005  and  a2<0), 
where  the  cloud  of  points  is  thick  and  dense,  04  takes  small  values 
meaning  a  weak  variation  of  the  pycnocline  depth.  a2  is  negative. 
Looking  at  the  profiles  projected  on  the  second  EOF  with  negative 
coefficients  on  Fig.  7,  the  corresponding  profiles  are  smooth  and  typical 
of  a  well  stratified  water  column.  This  regime  represents  a  “rest  regime” 
with  no  internal  wave  events.  Looking  now  at  the  “branches”  of  the 
“crescent”  (|aq|>- 0.005  and  a2>0),  one  can  see  how  both  coefficients 
increase  in  a  correlated  way.  This  regime  represents  typically  the 
internal  wave  dynamics:  when  the  pycnocline  is  depressed  by  the 
internal  solitary  waves  (represented  by  the  variation  of  ah),  the 
depression  of  the  interface  generates  a  horizontal  circulation  with 
opposite  horizontal  currents  below  and  above  the  interface  of  the 
pycnocline  thus  modifying  the  density  gradient  in  the  pycnocline 
(represented  by  a2).  These  baroclinic  currents  are  added  to  the  external 
barotropic  current.  Once  the  internal  wave  train  has  passed,  this 
baroclinic  structure  disappears  and  we  are  back  in  the  first  regime. 

3.  Modelling 

In  order  to  confirm  the  crescent-shape  scatter  diagram  obtained 
from  the  data,  several  numerical  simulations  were  made  using  the 
internal  wave  Lamb  model  (Lamb,  1994). 

3.1  The  2.5  dimensional  internal  wave  Lamb  model 

The  Lamb  (1994)  model  is  used  here  for  simulating  the  generation 
and  propagation  of  internal  solitary  waves  in  the  Strait  of  Messina.  It  is  a 


fully  nonlinear  and  nonhydrostatic  model.  It  is  based  on  the 
incompressible  Boussinesq  equations  on  a  rotating  /-plane.  In  the  y- 
direction,  the  velocity  V  is  included  but  its  derivatives  with  respect  to 
the  y-coordinate  are  neglected  (hence,  the  designation  “2.5  dimen¬ 
sional”  representation).  The  x  and  z  coordinates  define  a  vertical  plane 
with  the  x  coordinate  laid  out  along  bathymetry  contours  in  the  centre 
of  the  strait.  The  equations  of  the  model  are: 

Vt  +  V.W  -jVx  k  =  -  VP  -  kpg,  (4) 

pt  +  V.Vp  =  0,  (5) 

V.V  =  0,  (6) 

where  viscosity  is  neglected  for  internal  solitary  wave  generation 
and  propagation  (Warn  Varnas  et  al.,  2007),  V(u,v,w)  is  the  velocity 
vector,  V  is  the  three-dimensional  vector  gradient  operator,  subscript 
t  denotes  the  time  derivative,  p  is  the  density,  P  is  the  pressure,  g  is 
the  gravitational  acceleration,  /  is  the  Coriolis  parameter  taken  as 
9.03x10“  5  s-1,  and  k  is  the  unit  vector  along  the  z-direction,  z  is 
downward.  In  the  three-dimensional  Eqs.  (4)— (6)  the  partial  deriva¬ 
tives  with  respect  to  y  are  neglected,  i.e.  d  /  dy(.)  =  0.  Thus,  Eqs.  (4)- 
(6)  are  equivalent  to  Eqs.  (la)  to  (Id)  in  Lamb  (1994).  Before  the 
equations  are  solved,  they  are  transformed  to  a  terrain  following 
coordinate  system  in  the  vertical,  which  leads  to  higher  vertical 
resolution  over  the  sill  of  the  strait.  The  equations  are  solved  on  a 
domain  bounded  below  by  the  topography  at  z  =  h(z)  and  above  by  a 
free  slip  lid  at  z  =  H,  H  being  the  deep-water  depth.  At  the  upper 
boundary  there  is  no  normal  flow  and  the  normal  derivative  of  the 
tangential  velocity  is  set  to  zero.  At  the  lower  boundary  rigid  conditions 
are  used,  and  all  velocities  are  set  to  zero.  The  right  boundary  (northern 
basin)  is  treated  as  an  outflow.  On  the  left  boundary  (southern  basin), 
the  flow  ifib  is  forced  by  specifying  the  semidiurnal  tidal  velocity  as 
Uib  =  Vr1sin(fi)t  +  6)  where  VT  is  the  semidiurnal  tidal  magnitude,  co  the 
tidal  frequency  and  8  the  phase  factor.  Reflection  occurs  on  the  left  and 
the  simulation  length  time  must  be  adequately  chosen  in  order  to  avoid 
the  perturbation  of  the  physical  processes  by  the  reflected  waves. 
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The  model  was  initialized  with  150  vertical  levels,  which  gives  a 
vertical  resolution  of  50  cm  at  the  sill  and  8  m  at  the  southern 
boundary.  The  along-flow  coordinate  ranged  from  a  100  km  south  to 
100  km  north  of  the  sill  of  the  strait.  Only  the  M2  tidal  component  is 
considered  for  initial  tidal  forcing,  as  this  is  the  main  tidal  component 
in  the  Strait  (Alpers  and  Salusti,  1983;  Brandt  et  al.,  1998).  Several 
simulations  were  made  by  varying  the  semidiurnal  tidal  magnitude 
Vt-  The  topography  is  represented  by  cos h(x,z)  functions  (Fig.  10).  The 
density  field  for  model  initialization  was  derived  from  the  thermistor 
chain  measurements  for  the  southern  boundary  and  CTD  surveys  for 
the  northern  boundary.  The  model  predictions  are  for  the  median  of 
the  strait.  Five  hypothetical  moorings  were  implemented  in  the  model 
(Fig.  10).  We  concentrated  our  analysis  on  “mooring  2”  in  the  south  of 
the  strait  as  it  corresponds  to  the  position  of  the  thermistor  chain  in 
the  COACH06  experiment.  “Mooring  2”  was  defined  with  10  instru¬ 
ments  ranging  from  15  m  depth  to  132  m  depth. 

3.2.  EOF  analysis  of  the  model  outputs 

As  was  done  with  the  thermistor  data,  an  EOF  analysis  was  applied 
to  the  numerical  density  outputs  at  “mooring  2”  to  analyze  internal 
waves.  We  stress  here  that  we  are  adopting  a  statistical  approach,  not 
at  all  a  deterministic  approach  focusing  on  time  of  release  of  the  bore 
or  internal  wave  trains  precise  occurrence.  Variations  induced  by 
internal  solitary  waves  are  tracked  by  EOF  modes.  The  numerical  runs 
were  made  in  order  to  validate  the  characterization  of  the  internal 
waves  made  in  paragraph  2.3. 

Several  runs  were  made  by  changing  the  initial  tidal  M2 
component  speed  (respectively  0.01,  0.02,  0.03,  0.04,  0.05,  0.09  and 
0.12  ms-1),  all  initialized  by  the  same  density  profile  extracted  from 
the  thermistor  data.  Here  the  amplitude  of  the  barotropic  tidal 
component  is  varied  to  match  the  propagation  speed  of  the  internal 
solitary  wave  trains. 

For  an  initial  tidal  speed  of  0.01  m  s-  \  the  pycnocline  is  depressed 
by  the  initial  barotropic  tide  but  hardly  any  bore  is  generated.  For  an 
initial  tidal  speed  of  0.09  m  s-1,  the  amplitude  of  the  waves  is 
unrealistic  (over  200  m!)  and  they  propagate  too  fast  (Fig.  11). 

Given  the  SAR  figures,  the  run  closer  to  the  observed  internal 
waves  is  the  one  with  an  initial  barotropic  tidal  speed  of  0.03  m  s-1. 


The  propagation  speed  of  the  wave  train  is  1.06  m  s-  \  which  is  close 
to  the  propagation  speed  of  1.00  m  s- 1  obtained  from  the  SAR  data.  At 
25-30  km  south  of  the  sill,  the  modelled  wave  train  evolves  as  2  to  3 
internal  solitary  waves  (Fig.  11),  which  corresponds  well  to  the 
number  of  waves  observed  at  the  same  distance  on  the  SAR  images 
(Fig.  3).  The  amplitude  of  the  waves  (50  m)  at  30  km  from  the  sill  is 
also  close  to  that  of  the  thermistor  data  (Fig.  5).  The  pycnocline 
displacements  exhibit  a  pronounced  signature  of  internal  solitary 
wave  trains.  The  mid  pycnocline  location  can  simulate  the  interface 
deviation  for  a  two-layer  model.  Close  to  the  sea  surface  there  is  a 
signature  of  internal  solitary  waves  in  phase  with  the  pycnocline 
signature.  Time  evolution  of  pycnocline  and  sea  surface  signatures  can 
both  be  compared  to  SAR  observations.  The  first  two  EOF  modes  are  also 
well  illustrated.  The  first  EOF  mode  is  predominant  when  looking  at  the 
propagating  internal  wave  train  between  30  and  40  km  away  from  the 
sill  where  all  isopycnals  behave  similarly  in  the  pycnocline.  The  second 
EOF  mode,  where  isopycnals  are  moving  in  opposite  directions  above 
and  below  the  pycnocline,  is  present  close  to  the  sill  and  up  to  15  km 
southwards.  It  is  important  to  notice  that  an  unrealistic  barotropic  tidal 
speed  (such  as  runs  with  0.01  m  s- 1  and  0.09  m  s- 1  tidal  speed)  leads  to 
wrong  internal  baroclinic  wave  features. 

Looking  now  at  the  first  two  EOFs  for  the  run  with  0.03  m  s-1 
barotropic  tidal  velocity,  it  is  seen  that  the  EOF  extrema  are  deeper  for 
the  model  (Fig.  6)  than  for  the  data.  At  the  southern  boundary,  the 
model  is  initialized  by  a  density  profile  measured  when  the  pycnocline 
is  deeper  at  the  beginning  of  the  thermistor  string  measurement 
period.  The  thermistor  string  EOFs  computed  for  the  whole  period  of 
the  measurement  (Fig.  5)  thus  lead  to  shallower  maxima  because  of 
the  pycnocline  rise.  Nonetheless,  the  data  EOFs  and  the  model  EOFs 
agree  well  in  amplitude  and  shape. 

4.  Discussion  of  the  simulation  results  and  comparison  with  the 
data  EOF  distribution 

4.1.  Distribution  of  model  expansion  coefficients 

To  study  the  consistency  of  data  expansion  coefficients  in  a  scatter 
diagram,  the  previous  methodology  was  applied  to  three  numerical 
simulations  (Fig.  12)  at  “mooring  2”:  “No  IW  run”  (initial  tidal  speed 


Topography  used  for  Lamb  model  and  positions  of  the  "moorings" 


South  (Ionian  Sea)  <-  Along  Strait  distance  (km)  ->  (Tyrrhenian  Sea)  North 

Fig.  10.  Topography  (black  thick  line)  along  the  strait  used  in  the  2.5D  Lamb  (1994)  model.  Five  moorings,  called  ‘Moo’  on  the  graph  (grey  diamond  lines),  were  implemented  in  the 
model.  Moorings  1, 2, 3,  5  have  10  instruments  on  the  vertical,  Mooring  4  has  only  five  instruments  as  it  was  implemented  at  the  sill  of  the  strait.  Moorings  2  and  5  correspond  to  the 
positions  of  the  thermistor  strings  deployed  during  the  COACH06  cruise  (north  basin  was  not  considered  in  this  study). 
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Density  field  -  No  IW  run  (0.01  m/s) 
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Fig.  11. 2D  density  fields  after  20  h  run  for  an  initial  M2  tidal  speed  of  0.01  m  s_  1  (top),  0.03  m  s_  1  (middle)  and  0.09  m  s_  1  (bottom).  Density  units  are  “sigma-t”.  The  area  considered  is 
the  one  illustrated  in  the  zoom  of  Fig.  10.  The  amplitude  of  the  internal  waves  is  highly  correlated  to  the  initial  barotropic  tidal  speed.  If  the  barotropic  tidal  speed  is  too  high,  the  internal 
wave  features  become  unrealistic  (bottom)  for  the  strait  of  Messina.  This  parameter  must  be  tuned  very  carefully  in  order  to  get  proper  internal  solitary  wave  events. 


0.01  ms-1),  “IW-run”  (initial  tidal  speed  0.03  ms-1)  and  “unrealistic 
IW”  run  (initial  tidal  speed  0.09  m  s- 1)  where  “IW”  stands  for  internal 
wave. 

In  the  absence  of  internal  waves  (top  in  Fig.  12),  the  cloud  of  points  is 
small  (black  dots)  and  centered  at  the  cloud  of  points  of  the  data  (grey 
dots).  It  has  a  weak  expansion  of  a j  and  very  little  variations  of  a2.  Thus 
the  main  displacement  of  the  pycnocline  is  a  limited  vertical  translation 
with  no  internal  wave  propagation.  This  is  confirmed  by  the  2D  density 
field  on  Fig.  11  (top).  The  depression  of  the  pycnocline  due  to  the  bore, 
close  to  the  sill,  is  still  observed  though  this  bore  does  not  disintegrate 
into  a  proper  internal  wave  train.  Small  vertical  oscillations  of  the 
pycnocline  are  still  generated  further  down  the  strait.  These  oscillations 
are  represented  by  the  variation  of  ol\{  —  0.0055 <a^  <0.004). 

As  the  prescribed  barotropic  tidal  speed  increases  to  0.3  m  s-  \  we 
observe  a  great  expansion  of  the  “crescent”.  The  increase  in  the  values  of 
the  expansion  coefficients  corresponds  to  the  presence  of  bores  and  their 
subsequent  disintegration  into  solitary  wave  trains.  The  middle  graph  in 
Fig.  12  shows  the  best  fit  case  for  a  barotropic  initial  tidal  forcing  of 


0.03  ms-1.  The  model  distribution  in  this  case  fits  well  the  data  cloud  of 
points.  This  is  confirmed  by  the  nice  internal  solitary  wave  train  observed 
on  Fig.  11  (middle).  With  an  initial  tidal  forcing  of  0.09  ms-1,  the 
expansion  coefficient  distribution  of  the  model  does  not  fit  that  of  the 
data  anymore.  The  first  mode  interval  tends  towards  a  wider  negative 
limit,  which  means  the  pycnocline  gets  too  deep.  The  second  mode  has 
also  wider  positive  variations  where  it  used  to  be  negative  in  the  centre 
left  of  the  crescent  shape  ( —  0.002  <aA  <  0),  meaning  that  the  pycnocline 
slope  keeps  varying  while  the  pycnocline  depth  is  above  its  mean 
position.  This  illustrates  a  longer  time  needed  to  return  to  equilibrium  of 
the  pycnocline  as  the  internal  wave  characteristics  are  unrealistically 
large  (propagation  speed  of  1.30  m  s-1  with  wave  amplitude  of  200  m 
on  Fig.  11  bottom)  compared  with  those  of  the  data  (propagation  speed 
of  1.00  m  s- 1  with  waves  amplitude  of  50  m). 

These  results  confirm  that  incorrect  physics  leads  to  an  incorrect 
expansion  coefficient  distribution,  and  thus  that  the  scatter  diagram  is 
a  means  to  identify  and  characterize  internal  solitary  wave  trains 
occurring  in  the  area. 
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Fig.  12.  Comparisons  of  the  scatter  diagrams  (04  and  a2  versus  one  another)  of  the  thermistor  string  and  the  model.  The  grey  dots  represent  the  scatter  diagram  obtained  with  the 
expansion  coefficients  of  the  density  from  the  thermistor  string.  The  black  dots  represent  the  scatter  diagram  from  the  model  density  outputs  for  the  ‘No  IW’  simulation  (top),  the  ‘IW’ 
simulation  (middle)  and  the  ‘Unrealistic  IW’  simulation  (bottom),  corresponding  to  the  three  simulations  displayed  in  Fig.  11. 


4.2.  Model  scatter  diagram  sensitivity  to  initial  density  profile 

To  study  the  sensitivity  of  the  model  we  reconstruct  initial  density 
profiles  using  the  EOF  projection  coefficients  of  Fig.  9,  taken  inside  and 
outside  the  cloud  of  points.  We  then  reconstruct  the  corresponding  density 
profiles  using  Eq.  (2).  These  profiles  are  called  profile  A  (cq  =  —0.0071, 
a2  =  0),  profile  B  (cq  =  0,  a2  =  0.0041 ),  profile  C  (04  =  0.0048,  a2  =  0)  and 
profile  D  (04  =  0,  a2  =  —  0.0018).  They  are  shown  in  Fig.  9  by  their  EOF 
coefficients  and  are  outlined  in  Fig.  7.  As  these  profiles  are  not  reconstructed 
by  the  full  combination  of  all  EOF  modes,  they  are  not  necessarily  realistic, 
as  the  strong  density  inversion  in  profile  B  shows,  but  this  will  give  hints  on 
the  model  behaviour  and  stability. 

Profiles  A  and  C  are  reconstructed  with  an  anomaly  perturbation 
that  represents  an  up  and  down  movement  of  the  pycnocline  (Fig.  7 
left).  On  profile  A,  the  pycnocline  (inflexion  point)  is  higher  (70  m) 
than  the  mean  pycnocline  depth  (95  m).  On  the  contrary,  profile  C, 
with  a  high  positive  a !  coefficient  has  a  pycnocline  depth  deeper  than 
the  mean  profile.  For  these  two  profiles,  a2  is  null,  meaning  that  we 
have  filtered  out  the  second  mode. 

For  profiles  B  and  D,  04  is  null.  The  profiles  are  reconstructed  using 
only  the  second  EOF,  meaning  that  the  pycnocline  depth  corresponds 


to  the  mean  value  observed  in  the  data  and  that  the  pycnocline 
vertical  translation  is  filtered  out  (no  EOF1 ).  For  profile  B,  an  unstable 
slope  perturbation  of  the  pycnocline  is  introduced:  a2  has  a  very  high 
positive  value,  characteristic  of  a  strong  steepening  of  the  slope  of  the 
density  gradient  leading  to  density  inversion  (Fig.  7).  EOF2  represents 
only  4.5%  of  the  variability  as  mentioned  in  Section  3.2.  Reconstructing 
the  profiles  by  filtering  the  first  or  second  EOF,  specially  the  first  EOF 
(95%  of  the  variability),  may  not  lead  to  stable  profiles,  which  is  the 
case  for  this  latter  case. 

Profiles  A,  B  and  C  are  taken  out  of  the  cloud  of  points  (Fig.  9).  The 
EOF  coefficients  of  profile  D  have  been  chosen  in  order  to  locate 
the  density  profile  in  the  bottom  and  middle  of  the  cloud  of  points. 
The  reconstructed  profile  is  characteristic  of  a  stable  stratified  water 
column  with  a  mean  pycnocline  depth  at  95  m  and  a  strong  gradient 
(small  pycnocline  thickness).  The  aim  is  to  study  the  reaction  of  the 
model  to  initialisation  from  density  profiles  that  are  different  from 
the  measurements.  The  scatter  diagrams  for  each  profile  are  plotted  in 
Fig.  13. 

Our  first  remark  is  that  wherever  the  initial  profile  has  been  chosen 
inside  or  outside  the  data  distribution  of  points,  the  expansion  of  the 
model  scatter  diagram  fits  into  that  of  the  data. 
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Fig.  13.  Comparisons  of  data/modei  scatter  diagrams.  The  grey  dots  represent  the  scatter  diagram  from  the  thermistor  string  data.  The  black  dots  represent  the  scatter  diagram  from 
the  model  outputs  for  a  run  initialized  by  density  profile  A  (top  left),  by  density  profile  C  (top  right),  by  density  profile  B  (bottom  right)  and  by  density  profile  D  (the  four  profiles  are 
highlighted  on  Fig.  9).  A,  B  and  C  are  outside  the  cloud  of  points.  The  best  data/model  fit  is  for  the  run  initialized  by  profile  D  which  is  close  to  the  bottom  of  the  “crescent”. 


Profile  B  is  representative  of  a  density  profile  with  a  very  strong 
shear  leading  to  a  density  inversion.  The  corresponding  cloud  of 
points  looks  random  with  no  specific  shape  (Fig.  13  bottom  left).  The 
distribution  is  centred  around  the  origin  with  excursions  of  that 
reflect  small  depth  variations:  the  density  inversion  leads  to  instability 
and  creates  small  oscillations  of  the  pycnocline  interface 
(  —  0.004  <a^<  0.003)  with  a  very  weak  shear  (  —  0.001  <a2<  0.001). 
No  crescent  type  of  curvature  is  present  in  this  case  showing  that  the 
modeled  density  evolves  around  the  mean  density  profile.  The  initial 
density  inversion  is  smoothed  by  model  dynamics.  The  4.5% 
perturbation  input  through  EOF2  does  not  lead  to  a  random  behaviour 
of  the  model.  The  model  is  very  stable  and  sticks  to  the  average  profile. 
The  weak  shear  means  an  absence  of  internal  solitary  wave  events. 

Initializing  the  model  with  a  profile  whose  pycnocline  is  less  deep, 
profile  A,  results  in  the  scatter  diagram  shown  in  Fig.  13  top  left.  The 
variation  consists  of  line  movements  with  changing  shapes  and 
direction.  A  light  crescent  shape  starts  appearing.  There  is  still  some 
random  behavior  around  the  mean  position.  For  profiles  C  and  D,  the 
distribution  of  the  coefficients  contains  crescent  type  curvatures 
(Fig.  13  right,  top  and  bottom),  similar  to  the  variations  of  the  tuned 
model  distributions  in  Fig.  12.  The  gradient  around  the  inflexion  point 
of  profile  D  is  stronger  than  the  gradient  of  profile  C  (Fig.  7),  meaning 
that  the  pycnocline  is  thicker  for  profile  C  (about  100  m)  than  for 
profile  D  (about  70  m).  It  is  observed  on  the  scatter  diagram  of  profile 
C  that  oil  is  not  reaching  the  extreme  values  of  the  data  coefficient 
value  interval  although  the  initial  downward  displacement  of  the 
pycnocline  for  profile  C  results  in  a  scatter  diagram  that  is  in  range  of 
observed  scatter  diagram.  There  is  an  internal  wave  train,  but  as  the 
pycnocline  is  thicker,  the  amplitude  of  the  vertical  displacement  is  not 
as  strong  as  for  profile  D  (5  to  10  m  difference,  not  shown).  Another 
consequence  of  the  gradient  difference  is  that,  the  slope  being 
smoother  for  profile  C,  the  horizontal  circulation  generated  by  the 


interface  slope  will  be  weaker  and  thus  the  shear  will  not  be  as  strong: 
the  internal  solitary  waves  have  a  steeper  front  for  profile  D.  The 
smoother  pycnocline  gradient  of  profile  C  leads  to  “weaker”  internal 
waves  (smaller  amplitude)  than  the  observed  ones.  Profile  D  is  the 
profile  that  generates  the  internal  wave  train  closer  to  the  observed 
characteristics.  It  is  also  the  one  whose  distribution  fits  best  the  data 
scatter  diagram. 

All  four  of  the  constructed  initial  profiles,  Fig.  9,  lead  to  scatter 
diagrams  that  are  in  range  of  observations  (Fig.  13).  This  suggests  that 
variations  of  initial  pycnocline  depth  and  slope  of  pycnocline  are 
tolerable  in  model  initialization.  The  agreement  of  the  model 
predicted  scatter  diagrams  with  the  observations  shows  that  the 
Lamb  (1994)  model  represents  well  the  physics  and  phenomena 
taking  place. 

5.  Summary  and  conclusions 

The  Strait  of  Messina  with  its  shallow  sill  and  strong  tidal  currents 
leading  to  strong  interfacial  depressions  is  an  important  place  of 
internal  bore  generation  and  internal  solitary  wave  occurrence.  The 
nonlinear,  non-hydrostatic  Lamb  model,  initialized  by  two  density 
profiles  on  each  side  of  the  strait,  behaves  very  well  while  reproducing 
internal  wave  characteristics  and  is  in  a  very  good  qualitative  and 
quantitative  agreement  with  SAR  and  in-situ  data.  The  agreement  of 
model  with  measurements  indicates  that  the  Lamb  (1994)  model 
represents  well  the  physics  and  variability  of  the  internal  wave 
dynamics  and  phenomena  taking  place. 

In  this  study,  we  have  presented  an  original  method  for  internal 
wave  detection  in  the  south  of  the  Strait  of  Messina.  This  character¬ 
ization  is  based  on  EOF  analysis.  It  consists  in  studying  the  scatter 
diagram  formed  by  the  expansion  coefficient  of  the  first  EOF  mode 
versus  the  expansion  coefficient  of  the  second  EOF  mode.  When  there 


S204 


G.  Casagrande  et  al.  /  Journal  of  Marine  Systems  78  (2009)  S191-S204 


is  no  internal  wave,  the  scatter  diagram  is  small  and  the  projection  on 
the  first  mode  is  less  important  than  the  one  on  the  second  mode. 
When  there  is  an  internal  wave  event,  the  intervals  of  the  expansion 
coefficients  become  wider  so  that  the  scatter  diagram  assumes  a 
crescent  shape.  The  EOF  expansion  coefficients  and  their  scatter 
diagram  reflect  the  presence  of  internal  bores  and  internal  waves.  As 
the  scatter  diagram  has  a  well  defined  crescent  shape  distribution,  a 
polynomial  can  be  fitted  to  the  cloud  of  points  and  thus  give  a  relation 
between  the  two  expansion  coefficients.  The  validity  of  this  metho¬ 
dology  was  tested  by  comparing  model  outputs  with  in-situ  data.  The 
model  was  initialized  with  a  density  profile  derived  from  thermistor 
chain  and  CTD  probe  measurements.  The  tuning  was  accomplished 
by  varying  the  tidal  forcing.  The  predicted  internal  wave  speed  and 
amplitude  were  matched  with  SAR  and  in-situ  observations.  The 
scatter  diagram  of  the  model  output  and  the  scatter  diagram  derived 
from  measurements  are  in  good  agreement.  It  means  that  the  model 
can  predict  well  the  amplitude  of  the  two  main  EOF  modes,  but  also 
the  difference  in  phase  between  these  modes. 

This  methodology  could  be  interesting  for  operational  oceano¬ 
graphy  or  tactical  applications  in  navy  needs.  If  the  scatter  diagram  is 
well  shaped  as  observed  in  the  Strait  of  Messina,  a  polynomial 
function  can  be  fitted  to  the  crescent  shape  and  could  be  derived  for 
various  regions  of  interest.  This  polynomial  represents  the  relation 
between  the  first  two  expansion  coefficients.  It  could  be  sampled  in 
order  to  reconstruct  density  or  sound  speed  vertical  profiles  using  the 
first  two  EOF  vertical  modes.  The  fitted  function  is  then  a  way  to 
represent  all  the  water  column  profiles/physics  that  can  be  encoun¬ 
tered  in  the  area.  Once  this  function  is  established,  there  is  no  more 
need  for  in-situ  measurements  or  highly  demanding  computational 
cost  modelling.  Navy  forces  could  then  simulate  any  acoustical 
scenario  they  wish  to  study  for  any  chosen  profile,  or  may  combine 
the  curve  fitted  to  the  scatter  diagram  with  systematic  application  of 
acoustical  propagation  models  and  thus  get  a  complete  acoustical 
variability  map  of  the  area. 
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